import psycopg2
import psycopg2.extras
import sys
import gdal
from gdalconst import *
import os
import math
from shapely.wkt import dumps, loads
import datetime
from datetime import timedelta
import numpy as np
from constant import *
import constant


def interpolate(x1,y1,x0,y0,x):
    y=math.exp((x-x0)*((math.log(y1)-math.log(y0))/(x1-x0))+math.log(y0))
    return y

def corel_Npp_aeronet(in_file,time_range,neighbor_pixel):
    station_query = 'SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE id between 68 and 74'
    #mod_query="SELECT aqstime, st_neighborhood(rasref, (st_GeomFromText('POINT({0} {1})', 4326)), 4, 4, true) as aod FROM apom.satresampviirs as mod04 order by mod04.aqstime asc"
    mod_query = "SELECT aqstime, st_neighborhood(rasref, (st_GeomFromText('POINT({0} {1})', 4326)), {2}, {2}, true) as aod FROM apom.satresampviirs as viirs_aot order by viirs_aot.aqstime asc"
    aeronet_query = "SELECT avg(aod_675) as aod675,avg(aod_500) as aod500 from apom.grdaeraot_data where processtype = 2 and (aod_500>=0 and aod_675>=0) and aqstime > '{0}'::timestamp - interval '{1} minutes' and aqstime < '{0}'::timestamp + interval '{1} minutes' and stationid = {2} GROUP BY aqstime"
    
    f = open(in_file, "w+")
    con = None
    try:
        # doc danh sach tram
        con = psycopg2.connect(host=constant.__HOST__, database=constant.__DBNAME__, user=constant.__USERNAME__, password=constant.__PASSWORD__) 
        cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
       
        cur.execute(station_query)
        stations = cur.fetchall()
        for station in stations:
            stationid = station["id"]                
            location = loads(station["location"])
            # lay cac pixel xung quanh tram tren anh mod04,07
       
            cur.execute(mod_query.format(location.x, location.y,neighbor_pixel))
            rows = cur.fetchall()
            for row in rows:
                aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S")
                #aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
    
                # du lieu gia ve la 1 mang 2 chieu
                aod = row["aod"]
                aod_result = []
                for i in range(len(aod)):
                    # filter nhung gia tri != None thi se them vao mang ket qua
                    tmp = filter(lambda a: a != None, aod[i])
                    aod_result += tmp
                # calculate real aod value
                mod04_aod = (np.median(aod_result)) if len(aod_result) > 0 else 0
                    
                if mod04_aod > 0:
                    
                    cur.execute(aeronet_query.format(aqstime, time_range, stationid))                  
                    pm_row = cur.fetchone()
                    if not(pm_row is None):
                        aronet_aod=interpolate(675, pm_row["aod675"], 500, pm_row["aod500"], 550)
                        f.write("{0},{1},{2},{3},{4},{5}\n".format(stationid, aqstime.year, aqstime.month, aqstime, aronet_aod, mod04_aod))
                #print aqstime
    except psycopg2.DatabaseError, e:
        print 'Error %s' % e    
        sys.exit(1)
    finally:
        if con:
            con.close()
        if f:
            f.close()
    data = np.loadtxt(in_file, delimiter = ",", usecols=(4,5))
    a= np.corrcoef(data[:,0],data[:,1])
    print "number of sample",len(data)
    print "npp-aeronet R2:",a[1,0]

def corel_Npp_pm(in_file,time_range,neighbor_pixel):
    
    station_query='SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE id >=1 and id <=5'
    #mod_query="SELECT mod04.aqstime, mod04.aod from (SELECT aqstime, st_neighborhood(rasref, (st_GeomFromText('POINT({0} {1})', 4326)), 4, 4, true) as aod FROM apom.satresampviirs where aqstime < '2014-09-10 00:00:00'::timestamp) as mod04 order by mod04.aqstime asc"
    mod_query="SELECT aqstime, st_neighborhood(rasref, (st_GeomFromText('POINT({0} {1})', 4326)), {2}, {2}, true) as aod FROM apom.satresampviirs order by aqstime asc"
   
    pm_query="SELECT avg(pm1), avg(pm25), avg(pm10) FROM apom.grdpmhiscem_data WHERE (pm1 > 0 OR pm25 > 0 OR pm10 > 0) and aqstime > '{0}'::timestamp - interval '{1} minutes' and aqstime < '{0}'::timestamp + interval '{1} minutes' and stationid = {2} GROUP BY aqstime"
    
    f = open(in_file, "w+")
    con = None
    try:
        # doc danh sach tram
        con = psycopg2.connect(host=constant.__HOST__, database=constant.__DBNAME__, user=constant.__USERNAME__, password=constant.__PASSWORD__) 
        cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
        cur.execute(station_query)
        stations = cur.fetchall()
        for station in stations:
            stationid = station["id"]
            # bo qua tram Khanh Hoa
            # if stationid  == 3:
            #     continue
                
            location = loads(station["location"])
    
            cur.execute(mod_query.format(location.x, location.y,neighbor_pixel))
         
            rows = cur.fetchall()
            for row in rows:
                aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
            
                # du lieu gia ve la 1 mang 2 chieu
                aod = row["aod"]
                aod_result = []
                for i in range(len(aod)):
                    # filter nhung gia tri != None thi se them vao mang ket qua
                    tmp = filter(lambda a: a != None, aod[i])
                    aod_result += tmp
                # calculate real aod value
                aod_value = (np.median(aod_result)) if len(aod_result) > 0 else 0
                if (aod_value > 0):
                    cur.execute(pm_query.format(aqstime, time_range, stationid))
                    pm_row = cur.fetchone()
                    if not(pm_row is None):
                        pm1 = pm_row[0]
                        pm25 = pm_row[1]
                        pm10 = pm_row[2]                 
                        f.write("{0},{1},{2},{3},{4},{5},{6},{7}\n".format(stationid, aqstime.year, aqstime.month, aqstime, pm1, pm25, pm10, aod_value))
            
    except psycopg2.DatabaseError, e:
        print 'Error %s' % e    
        sys.exit(1)
    finally:
        if con:
            con.close()
        if f:
            f.close()
    data = np.loadtxt(in_file, delimiter = ",", usecols=(4,5,6,7))
    print "Number of sample:",len(data)
    
    a= np.corrcoef([data[:,0],data[:,3]])
    print "NPP - PM1 R2:", a[1,0]
     
    b= np.corrcoef([data[:,1],data[:,3]])
    print "NPP - PM25 R2:", b[1,0]
     
    b= np.corrcoef([data[:,2],data[:,3]])
    print "NPP - PM10 R2:", b[1,0]


def corel_NppTemp_pm(in_file,time_range,neighbor_pixel):
    
    station_query = 'SELECT id, ST_AsText(location) as location FROM public.grdstation WHERE id >=1 and id <=5'
    #mod_query="SELECT aqstime, ST_Value(rasref, (st_GeomFromText('POINT({0} {1})', 4326)),true) as temp FROM apom.satresampmod07temperature where aqstime < '2014-09-10 00:00:00'::timestamp"
    mod_query = "SELECT aqstime, st_neighborhood(rasref, (st_GeomFromText('POINT({0} {1})', 4326)), {2}, {2}, true) as temp FROM apom.satresampviirstemperature order by aqstime asc"
    
    pm_query = "SELECT avg(pm1), avg(pm25), avg(pm10),avg(temp) FROM apom.grdpmhiscem_data WHERE (pm1 > 0 OR pm25 > 0 OR pm10 > 0) and aqstime > '{0}'::timestamp - interval '{1} minutes' and aqstime < '{0}'::timestamp + interval '{1} minutes' and stationid = {2} GROUP BY aqstime"
   
    f = open(in_file, "w+")
    con = None
    try:
        # doc danh sach tram
        con = psycopg2.connect(host=constant.__HOST__, database=constant.__DBNAME__, user=constant.__USERNAME__, password=constant.__PASSWORD__) 
        cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
        cur.execute(station_query)
        stations = cur.fetchall()
        for station in stations:
            stationid = station["id"]
            # bo qua tram Khanh Hoa
            # if stationid  == 3:
            #     continue
                
            location = loads(station["location"])
            cur.execute(mod_query.format(location.x, location.y,neighbor_pixel))                                                                                                                                                                                                                  
            rows = cur.fetchall()
            
            for row in rows:
                temp_aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
                temp = row["temp"]
                
                
                temp_result = []
                for i in range(len(temp)):
                    # filter nhung gia tri != None thi se them vao mang ket qua
                    tmp = filter(lambda a: a != None, temp[i])
                    temp_result += tmp
                # calculate real aod value
                temp_value = (np.median(temp_result)) if len(temp_result) > 0 else 0
                
                """
                if(temp!=None):
                    temp_value = temp
                else:
                    temp_value=0;
                 """
               
                if temp_value > 0:
                 
                    cur.execute(pm_query.format(temp_aqstime, time_range, stationid))
                    
                    pm_row = cur.fetchone()
                    if not(pm_row is None):
                        pm1 = pm_row[0]
                        pm25 = pm_row[1]
                        pm10 = pm_row[2]
                        cem_temp=pm_row[3]
                      
                        f.write("{0},{1},{2},{3},{4},{5},{6},{7},{8}\n".format(stationid, temp_aqstime.year, temp_aqstime.month, temp_aqstime, pm1, pm25, pm10,cem_temp,temp_value))
            
    except psycopg2.DatabaseError, e:
        print 'Error %s' % e    
        sys.exit(1)
    finally:
        if con:
            con.close()
        if f:
            f.close()
    data = np.loadtxt(in_file, delimiter = ",", usecols=(4,5,6,7,8))
    print "Number of sample:",len(data)
    
    a= np.corrcoef([data[:,0],data[:,4]])
    print "temp - PM1 R2:", a[1,0]
     
    b= np.corrcoef([data[:,1],data[:,4]])
    print "temp - PM25 R2:", b[1,0]
     
    b= np.corrcoef([data[:,2],data[:,4]])
    print "temp - PM10 R2:", b[1,0]
    
    b= np.corrcoef([data[:,3],data[:,4]])
    print "temp - CEM temp R2:", b[1,0]

 
if __name__ == '__main__':
    
    corel_NppTemp_pm("correl/npptemp_pm_30s_6km.csv",30,1)
    corel_NppTemp_pm("correl/npptemp_pm_30s_12km.csv",30,2)
    corel_NppTemp_pm("correl/npptemp_pm_30s_18km.csv",30,3)
    corel_NppTemp_pm("correl/npptemp_pm_30s_24km.csv",30,4)
    corel_NppTemp_pm("correl/npptemp_pm_30s_30km.csv",30,5)
    corel_NppTemp_pm("correl/npptemp_pm_30s_36km.csv",30,6)
    corel_NppTemp_pm("correl/npptemp_pm_30s_42km.csv",30,7)
    corel_NppTemp_pm("correl/npptemp_pm_30s_48km.csv",30,8)
    corel_NppTemp_pm("correl/npptemp_pm_30s_54km.csv",30,9)
 
    
    corel_NppTemp_pm("correl/npptemp_pm_60s_6km.csv",60,1)
    corel_NppTemp_pm("correl/npptemp_pm_60s_12km.csv",60,2)
    corel_NppTemp_pm("correl/npptemp_pm_60s_18km.csv",60,3)
    corel_NppTemp_pm("correl/npptemp_pm_60s_24km.csv",60,4)
    corel_NppTemp_pm("correl/npptemp_pm_60s_30km.csv",60,5)
    corel_NppTemp_pm("correl/npptemp_pm_60s_36km.csv",60,6)
    corel_NppTemp_pm("correl/npptemp_pm_60s_42km.csv",60,7)
    corel_NppTemp_pm("correl/npptemp_pm_60s_48km.csv",60,8)
    corel_NppTemp_pm("correl/npptemp_pm_60s_54km.csv",60,9)
    
    corel_NppTemp_pm("correl/npptemp_pm_90s_6km.csv",90,1)
    corel_NppTemp_pm("correl/npptemp_pm_90s_12km.csv",90,2)
    corel_NppTemp_pm("correl/npptemp_pm_90s_18km.csv",90,3)
    corel_NppTemp_pm("correl/npptemp_pm_90s_24km.csv",90,4)
    corel_NppTemp_pm("correl/npptemp_pm_90s_30km.csv",90,5)
    corel_NppTemp_pm("correl/npptemp_pm_90s_36km.csv",90,6)
    corel_NppTemp_pm("correl/npptemp_pm_90s_42km.csv",90,7)
    corel_NppTemp_pm("correl/npptemp_pm_90s_48km.csv",90,8)
    corel_NppTemp_pm("correl/npptemp_pm_90s_54km.csv",90,9)
       
    corel_NppTemp_pm("correl/npptemp_pm_720s_6km.csv",720,1)
    corel_NppTemp_pm("correl/npptemp_pm_720s_12km.csv",720,2)
    corel_NppTemp_pm("correl/npptemp_pm_720s_18km.csv",720,3)
    corel_NppTemp_pm("correl/npptemp_pm_720s_24km.csv",720,4)
    corel_NppTemp_pm("correl/npptemp_pm_720s_30km.csv",720,5)
    corel_NppTemp_pm("correl/npptemp_pm_720s_36km.csv",720,6)
    corel_NppTemp_pm("correl/npptemp_pm_720s_42km.csv",720,7)
    corel_NppTemp_pm("correl/npptemp_pm_720s_48km.csv",720,8)
    corel_NppTemp_pm("correl/npptemp_pm_720s_54km.csv",720,9)
    
    
    """
    corel_Npp_pm("correl/npp_pm_30s_6km.csv",30,1)
    corel_Npp_pm("correl/npp_pm_30s_12km.csv",30,2)
    corel_Npp_pm("correl/npp_pm_30s_18km.csv",30,3)
    corel_Npp_pm("correl/npp_pm_30s_24km.csv",30,4)
    corel_Npp_pm("correl/npp_pm_30s_30km.csv",30,5)
    corel_Npp_pm("correl/npp_pm_30s_36km.csv",30,6)
    corel_Npp_pm("correl/npp_pm_30s_42km.csv",30,7)
    corel_Npp_pm("correl/npp_pm_30s_48km.csv",30,8)
    corel_Npp_pm("correl/npp_pm_30s_54km.csv",30,9)
 
    
    corel_Npp_pm("correl/npp_pm_60s_6km.csv",60,1)
    corel_Npp_pm("correl/npp_pm_60s_12km.csv",60,2)
    corel_Npp_pm("correl/npp_pm_60s_18km.csv",60,3)
    corel_Npp_pm("correl/npp_pm_60s_24km.csv",60,4)
    corel_Npp_pm("correl/npp_pm_60s_30km.csv",60,5)
    corel_Npp_pm("correl/npp_pm_60s_36km.csv",60,6)
    corel_Npp_pm("correl/npp_pm_60s_42km.csv",60,7)
    corel_Npp_pm("correl/npp_pm_60s_48km.csv",60,8)
    corel_Npp_pm("correl/npp_pm_60s_54km.csv",60,9)
    
    corel_Npp_pm("correl/npp_pm_90s_6km.csv",90,1)
    corel_Npp_pm("correl/npp_pm_90s_12km.csv",90,2)
    corel_Npp_pm("correl/npp_pm_90s_18km.csv",90,3)
    corel_Npp_pm("correl/npp_pm_90s_24km.csv",90,4)
    corel_Npp_pm("correl/npp_pm_90s_30km.csv",90,5)
    corel_Npp_pm("correl/npp_pm_90s_36km.csv",90,6)
    corel_Npp_pm("correl/npp_pm_90s_42km.csv",90,7)
    corel_Npp_pm("correl/npp_pm_90s_48km.csv",90,8)
    corel_Npp_pm("correl/npp_pm_90s_54km.csv",90,9)
       
    corel_Npp_pm("correl/npp_pm_720s_6km.csv",720,1)
    corel_Npp_pm("correl/npp_pm_720s_12km.csv",720,2)
    corel_Npp_pm("correl/npp_pm_720s_18km.csv",720,3)
    corel_Npp_pm("correl/npp_pm_720s_24km.csv",720,4)
    corel_Npp_pm("correl/npp_pm_720s_30km.csv",720,5)
    corel_Npp_pm("correl/npp_pm_720s_36km.csv",720,6)
    corel_Npp_pm("correl/npp_pm_720s_42km.csv",720,7)
    corel_Npp_pm("correl/npp_pm_720s_48km.csv",720,8)
    corel_Npp_pm("correl/npp_pm_720s_54km.csv",720,9)
    
    corel_Npp_aeronet("correl/npp_aronet_60s_6km.csv",60,1)
    corel_Npp_aeronet("correl/npp_aronet_60s_12km.csv",60,2)
    corel_Npp_aeronet("correl/npp_aronet_60s_18km.csv",60,3)
    corel_Npp_aeronet("correl/npp_aronet_60s_24km.csv",60,4)
    corel_Npp_aeronet("correl/npp_aronet_60s_30km.csv",60,5)
    corel_Npp_aeronet("correl/npp_aronet_60s_36km.csv",60,6)
    corel_Npp_aeronet("correl/npp_aronet_60s_42km.csv",60,7)
    corel_Npp_aeronet("correl/npp_aronet_60s_48km.csv",60,8)
    corel_Npp_aeronet("correl/npp_aronet_60s_54km.csv",60,9)
    
    corel_Npp_aeronet("correl/npp_aronet_90s_6km.csv",90,1)
    corel_Npp_aeronet("correl/npp_aronet_90s_12km.csv",90,2)
    corel_Npp_aeronet("correl/npp_aronet_90s_18km.csv",90,3)
    corel_Npp_aeronet("correl/npp_aronet_90s_24km.csv",90,4)
    corel_Npp_aeronet("correl/npp_aronet_90s_30km.csv",90,5)
    corel_Npp_aeronet("correl/npp_aronet_90s_36km.csv",90,6)
    corel_Npp_aeronet("correl/npp_aronet_90s_42km.csv",90,7)
    corel_Npp_aeronet("correl/npp_aronet_90s_48km.csv",90,8)
    corel_Npp_aeronet("correl/npp_aronet_90s_54km.csv",90,9)
    
    corel_Npp_aeronet("correl/npp_aronet_720s_6km.csv",720,1)
    corel_Npp_aeronet("correl/npp_aronet_720s_12km.csv",720,2)
    corel_Npp_aeronet("correl/npp_aronet_720s_18km.csv",720,3)
    corel_Npp_aeronet("correl/npp_aronet_720s_24km.csv",720,4)
    corel_Npp_aeronet("correl/npp_aronet_720s_30km.csv",720,5)
    corel_Npp_aeronet("correl/npp_aronet_720s_36km.csv",720,6)
    corel_Npp_aeronet("correl/npp_aronet_720s_42km.csv",720,7)
    corel_Npp_aeronet("correl/npp_aronet_720s_48km.csv",720,8)
    corel_Npp_aeronet("correl/npp_aronet_720s_54km.csv",720,9)
    
    """

   